package scu.maqiang.cte;

import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3H201;
import scu.maqiang.fes.FES3H202;
import scu.maqiang.mesh.Mesh3H20;
import scu.maqiang.mesh.Mesh3H8;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.mesh.Tecplot;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class BrickResultH20_SC2 {

	public static void main(String[] args) {
		double kxx = 1.0, density = 1.0, specficHeat = 1.0;
		double E = 2.0, Nu = 0.3, alpha = 0.05;
		double beta = E / (1 - 2 * Nu) * alpha;

		double theta = 0.5, endTime = 10, dt = endTime / 1000;
		double omega = 0.25;

	    Mesh3H8 mesh3D0 = new Mesh3H8().cube(2, 2, 4);
	    mesh3D0.scale(0.5, 0.5, 1.0);

	    Mesh3H20 mesh3D = new Mesh3H20(mesh3D0);

		// �����¶ȳ���Ԫ�ռ�, ��װ��Ӧ�������
		FES3H201 fs = new FES3H201(mesh3D);
		SRMatrix LMatrix = new SRMatrix(fs.GetNdof());
		SRMatrix RMatrix = new SRMatrix(fs.GetNdof());
		ScalarFunc kxxFunc = (xy, llabel, param) -> kxx;
		ScalarFunc rhoFunc = (xy, llabel, param) -> density * specficHeat;
		fs.assembleHeatThetaScheme(new ScalarFunc[] { kxxFunc, rhoFunc }, new double[] { dt, theta },
				new BVPType[] { BVPType.COMMON, BVPType.CONSIST_MASS}, LMatrix, RMatrix);
		double[] x = new double[fs.GetNdof()];
		double[] RHS = new double[fs.GetNdof()];

		// ���õ��Գ�����Ԫ�ռ�, ��װ��Ӧ�������
		FES3H202 fs2 = new FES3H202(mesh3D);
		SRMatrix K = new SRMatrix(fs2.GetNdof());
		SRMatrix M = new SRMatrix(fs2.GetNdof());
		double[] d = new double[fs2.GetNdof()];
		double[] d0 = new double[fs2.GetNdof()];
		double[] vel = new double[fs2.GetNdof()];
		double[] vel0 = new double[fs2.GetNdof()];
		double[] acc = new double[fs2.GetNdof()];
		double[] acc0 = new double[fs2.GetNdof()];
		double[] heatLoad = new double[fs2.GetNdof()];
		double[] totalLoad = new double[fs2.GetNdof()];

		fs2.assembleStiff(new double[] { E, Nu }, BVPType.COMMON, K);
		fs2.assembleMass(new double[] { density }, BVPType.CONSIST_MASS, M);
		ScalarFunc myFunc = (coord, llabel, param) -> coord[2] > 0.99999 ? 1.0 : 0.0;
		double[] x0 = fs.valueFromFunc(myFunc, null);
		ScalarFunc betaFunc = (xy, llabel, param) -> beta;

		double[][] vecX = new double[1][];
		vecX[0] = x0;
		fs2.assembleHeatLoad(new ScalarFunc[] { betaFunc }, null, vecX, BVPType.COMMON, heatLoad);

		// 求解初始加速度
		NewIterSSolver solverM = new NewIterSSolver(M);
		solverM.PCGSSOR(heatLoad, acc, 1.5, 1);

		for (int i = 0, n = fs2.GetNdof(); i < n; i++) {
			acc0[i] = acc[i];
		}

		double[][] uvw = new double[3][fs.GetNdof()];
		fs2.extractComponent(acc0, uvw);

		String dir = "BrickResultH20\\";

		mesh3D.toTecplot(dir + "acc.dat", uvw);
		MVO.fill(uvw, 0.0);

		// K.sort();
		// System.out.println(K);
		// System.exit(0);

		SRMatrix ELMatrix = K.add((1.0 / (0.25 * dt * dt)), M);
		// ELMatrix.sort();
		// ELMatrix.toString();
		// System.out.println(ELMatrix);
		// System.exit(0);

		double time = 0.0;

		String heatoutputFile = dir + "Heat100_2.dat";
		mesh3D.toTecplot(heatoutputFile, time, x);

		String elasticOutputFile = dir + "Deform100_2.dat";
		mesh3D.toTecplotDeform(elasticOutputFile, 0.0, uvw, 0.0);
		// call mesh3D.exportMeshDeformTime(elasticOutputFile, time, x, x, 0.0D0)

		int lineNum = 100;
		double dis = 1.0 / lineNum;
		double[] lineX = new double[lineNum + 1];
		double[] lineY = new double[lineNum + 1];
		Arrays.setAll(lineX, i -> i * dis);

		String XYLineFile = dir + "Line100_2.dat";
		Tecplot.lineXYTime(XYLineFile, lineX, lineY, time);

		// x = 1.0D0
		String expansion = dir + "Expansion100_2.dat";
		String temp00File = dir + "temp00_100_2.dat";
		FileWriter fw = null;
		FileWriter fw2 = null;
		try {
			fw = new FileWriter(expansion);
			fw.write("0.0\t0.0\n");
			fw2 = new FileWriter(temp00File);
			fw2.write("0.0\t0.0\n");
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
		// XYStressYFile = trim(dir)//'StressYLine100_2.dat'
		// call exportXYTime(XYStressYFile, lineX, time, lineY)

		fs.applyBC_MBN(LMatrix, Direct.All, 2);

		fs2.applyBC_MBN(K, Direct.X, 3, 4, 5, 6);
		fs2.applyBC_MBN(K, Direct.Y, 3, 4, 5, 6);
		fs2.applyBC_MBN(K, Direct.Z, 1);

		fs2.applyBC_MBN(ELMatrix, Direct.X, 3, 4, 5, 6);
		fs2.applyBC_MBN(ELMatrix, Direct.Y, 3, 4, 5, 6);
		fs2.applyBC_MBN(ELMatrix, Direct.Z, 1);

		NewIterSSolver SolverLM = new NewIterSSolver(LMatrix);
		double[] Md = new double[fs2.GetNdof()];
		double[] MMd = Md.clone();

		NewIterSSolver solverEL = new NewIterSSolver(ELMatrix);
		time = 0.0;
		double coef = (1 - Nu) / (alpha * (1 + Nu));
		while (time < endTime) {
			time += dt;
			System.out.println("time: " + time);

			fs.applyBC_MBN(RHS, Direct.All, 1.0, 2);
			SolverLM.PCGSSOR(RHS, x, 1.5, 1);
			mesh3D.toTecplot(heatoutputFile, time, x);
			// RHS = RMatrix * x
			RMatrix.mul(x, RHS);

			for (int i = 0; i < lineNum + 1; i++) {
				lineY[i] = fs.getValue(x, new double[] { 0.25, 0.25, lineX[i] });
			}
			Tecplot.lineXYTime(XYLineFile, lineX, lineY, time);

			Arrays.fill(heatLoad, 0.0);
			vecX[0] = x;
			fs2.assembleHeatLoad(new ScalarFunc[] { betaFunc }, null, vecX, BVPType.COMMON, heatLoad);
			// call fs2.assembleHeatLoad(beta, x, 0, HeatLoad)
			// M.mul(1.0 / (0.25 * dt * dt), d, Md);
			Arrays.setAll(Md,
					i -> 1.0 / (0.25 * dt * dt) * d[i] + (1.0 / (0.25 * dt)) * vel[i] + (0.5 / 0.25 - 1.0) * acc[i]);
			M.mul(Md, MMd);
			for (int ii = 0, nn = fs2.GetNdof(); ii < nn; ii++) {
				totalLoad[ii] = heatLoad[ii] + MMd[ii];
			}
			// System.out.println(MVO.toString(totalLoad));
			// System.out.println(MVO.max(totalLoad) + "\t" + MVO.min(totalLoad));

			fs2.applyBCZero(totalLoad, Direct.X, 3, 4, 5, 6);
			fs2.applyBCZero(totalLoad, Direct.Y, 3, 4, 5, 6);
			fs2.applyBCZero(totalLoad, Direct.Z, 1);

			solverEL.PCGSSOR(totalLoad, d, 1.5, 1);
			fs2.extractComponent(d, uvw);
			mesh3D.toTecplotDeform(elasticOutputFile, time, uvw, coef);

			// call fs2.computeStrain(u, v, ex, ey, exy)
			// call fs2.computeStress(ex, ey, exy, u, E, Nu, tx, ty, txy, u, 1)

			// do i = 1, lineNum + 1
			// call fs.valueXY(ty, 0.25D0, lineX(i), val)
			// lineY(i) = val
			// end do
			// call exportXYTime(XYStressYFile, lineX, time, lineY)

			Arrays.setAll(acc, i -> 1 / (omega * dt * dt) * (d[i] - d0[i]) - 1 / (omega * dt) * vel0[i]
					- (0.5 / omega - 1) * acc0[i]);
			Arrays.setAll(vel, i -> vel0[i] + ((1 - 0.5) * acc0[i] + 0.5 * acc[i]) * dt);
			Arrays.setAll(acc0, i -> acc[i]);
			Arrays.setAll(vel0, i -> vel[i]);
			Arrays.setAll(d0, i -> d[i]);

			try {
				fw.write(time + "\t " + fs.getValue(uvw[2], new double[] { 0.25, 0.25, 1.0 }) * coef + "\n");
				fw2.write(time + "\t " + fs.getValue(x, new double[] { 0.0, 0.0, 0.0 }) + "\n");
			} catch (Exception e) {
				e.printStackTrace();
			}
		}
		try {
			fw.flush();
			fw.close();
			fw2.flush();
			fw2.close();
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
	}

}
